library(GEOquery)
library(tidyverse)
library(tidySummarizedExperiment)
library(ggpubr)

# GSE103740 human 2019 ------------
gset <- getGEO("GSE103740", AnnotGPL = TRUE) |>
  pluck(1) |>
  makeSummarizedExperimentFromExpressionSet()

gset@colData |>
  as_tibble() |>
  head() |>
  DT::datatable()

meta_res <- gset@colData |>
  as_tibble() |>
  mutate(infected = str_remove(characteristics_ch1.8, 'hiv infection: '),
         donor = str_remove(characteristics_ch1, 'donor: '),
         treatment = str_remove(characteristics_ch1.7, 'treatment: '),
         ttt = str_remove(characteristics_ch1.6, 'treatment: '),
         stim = str_remove(characteristics_ch1.3, 'stimulation: ')) |>
  select(title, geo_accession, infected, treatment, donor, ttt, stim)

meta_res |>
  summarise(n(),.by = time2)

tidy_res <- meta_res |>
  mutate(infected = ifelse(str_detect(treatment, 'Yes'), 'yes', 'no'),
         vaccinated = ifelse(ttt == 'VACCINE' | str_detect(treatment, 'vaccine'), 'yes', 'no'),
         time = str_extract(title, 'p.+vacc')) |>
  select(infected, vaccinated, donor, geo_accession, stim, time)

tidy_res |>
  summarise(n(),.by = c(time, donor))

## array annotation ------
probe_gene <- gset |>
  as_tibble()

hmces_probe <- probe_gene |>
  filter(`Gene symbol` == 'HMCES') |>
  select(1,3)

mat_hm <- gset |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  filter(ID == 'ILMN_1815682') |>
  pivot_longer(where(is.numeric))

mat_hm |>
  dplyr::rename(geo_accession = name) |>
  left_join(tidy_res) |>
  filter(time == 'post-vacc' & stim != 'DMSO') |>
  ggplot(aes(infected, value, color = infected)) +
  stat_summary(geom = 'crossbar', fun = 'mean') +
  geom_jitter(height = 0, width = .3) +
  stat_compare_means(method = 't.test', color ='black') +
  scale_color_manual(values = c('blue','red')) +
  scale_x_discrete(label = c('healthy','infected')) +
  labs(x = 'Infected after vaccination', y = 'HMCES expression') +
  theme_pubr()

# GSE72624 macaca 2016 ------------
gset <- getGEO("GSE72624", AnnotGPL = TRUE)

gset <- gset[[1]]

gset <- gset |> makeSummarizedExperimentFromExpressionSet()

meta_res <- gset@colData |>
  as_tibble() |>
  mutate(resistance = str_remove(characteristics_ch1.4, 'number of siv challenge to infection: ')) |>
  select(title, geo_accession, resistance)

tidy_res <- meta_res |>
  mutate(donor = str_extract(title, '(?<=donor ).+') |>
           str_remove('_.+'),
         time_point = str_extract(title, '.+(?= blood)'),
         resistance = ifelse(resistance == 'neg', '11', resistance) |> as.numeric())

## array annotation ------
probe_gene <- gset@featureData@data |>
  as_tibble()

hmces_probe <- probe_gene |>
  filter(`Gene symbol` == 'HMCES') |>
  select(1,3)

mat_hm <- gset |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  right_join(hmces_probe) |>
  pivot_longer(where(is.numeric))

mat_hm |>
  ggplot(aes(ID, value)) +
  geom_boxplot()

mat_hm |>
  filter(str_detect(ID, '1815')) |>
  dplyr::rename(geo_accession = name) |>
  left_join(tidy_res) |>
  mutate(time_point = fct_relevel(time_point, 'pre-vaccination')) |>
  ggplot(aes(resistance, value)) +
  geom_point() +
  geom_smooth(method = 'glm') +
  stat_cor() +
  facet_wrap(~time_point) +
  labs(x = 'resistance to infection', y = 'HMCES expression') +
  theme_pubr()

# GSE108011 macaca 2018 ------------
gset <- getGEO("GSE108011", AnnotGPL = TRUE) |>
  pluck(1) |>
  makeSummarizedExperimentFromExpressionSet()

gset@colData |>
  as_tibble() |>
  DT::datatable()

meta_res <- gset@colData |>
  as_tibble() |>
  mutate(resistance = str_remove(characteristics_ch1.4, 'number of siv challenge to infection: '),
         donor = str_remove(characteristics_ch1, 'donor: ')) |>
  select(title, geo_accession, resistance, donor)

meta_res |>
  summarise(n(),.by = donor)

tidy_res <- meta_res |>
  mutate(time_point = str_extract(title, 'prevax|24h.+|1we.+|2wee.+'),
         resistance = ifelse(resistance == 'neg', '7', resistance) |> as.numeric())

tidy_res |>
  summarise(n(),.by = time_point)

## array annotation ------
probe_gene <- gset@featureData@data |>
  as_tibble()

hmces_probe <- probe_gene |>
  filter(`Gene symbol` == 'HMCES') |>
  select(1,3)

mat_hm <- gset |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  right_join(hmces_probe) |>
  pivot_longer(where(is.numeric))

mat_hm |>
  filter(str_detect(ID, '1815')) |>
  dplyr::rename(geo_accession = name) |>
  left_join(tidy_res) |>
  mutate(time_point = fct_relevel(time_point, 'prevax')) |>
  ggplot(aes(resistance, value)) +
  geom_point() +
  geom_smooth(method = 'glm') +
  stat_cor() +
  facet_wrap(~time_point) +
  labs(x = 'resistance to infection', y = 'HMCES expression') +
  theme_pubr()
